if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
;!p.multi=[0,3,2]
;!p.charsize=3.5
!p.multi=[0,3,1] & !p.charsize=2.0 & !x.margin=[3.75,0] & !y.margin=[3,0.5]

default,nn,5
ntimes = (size(u))[3]
good = (ntimes-nn)+indgen(nn)-1
print,'doing averages over',nn, ' xaver records ...'

ruu = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
rvv = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
rww = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)


uu=fltarr(l,m)
vv=fltarr(l,m)
ww=fltarr(l,m)
oox=fltarr(l,m)
ooy=fltarr(l,m)
ooz=fltarr(l,m)

omega=fltarr(l,m)

quu=fltarr(l,m)
qvv=fltarr(l,m)
qww=fltarr(l,m)
qwv=fltarr(l,m)
qwu=fltarr(l,m)
qvu=fltarr(l,m)
rsinth=fltarr(l,m)
rurms2=fltarr(l,m)
frmc=fltarr(l,m)
frrs=fltarr(l,m)
fthmc=fltarr(l,m)
fthrs=fltarr(l,m)

frt=fltarr(l,m)
ftht=fltarr(l,m)
fx=fltarr(l,m)
fy=fltarr(l,m)

urms2t = ruu+rvv+rww
urms2 = (total(urms2t[good]/float((size(good))[1])))

; reynolds stresses
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=rr*r[k]*sin(theta[j])
    uu[k,j] = total(u[j,k,good])/float((size(good))[1])
    vv[k,j] = total(v[j,k,good])/float((size(good))[1])
    ww[k,j] = total(w[j,k,good])/float((size(good))[1])
    oox[k,j] = total(ox[j,k,good])/float((size(good))[1])
    ooy[k,j] = total(oy[j,k,good])/float((size(good))[1])
    ooz[k,j] = total(oz[j,k,good])/float((size(good))[1])
    quu[k,j]=rho[j,k]*total(u2[j,k,good]- u[j,k,good]^2)/float((size(good))[1])
    qvv[k,j]=rho[j,k]*total(v2[j,k,good]- v[j,k,good]^2)/float((size(good))[1])
    qww[k,j]=rho[j,k]*total(w2[j,k,good]- w[j,k,good]^2)/float((size(good))[1])
    rurms2[k,j]=quu[k,j]+qvv[k,j]+qww[k,j]
;
    qwv[k,j]=total(rwv[j,k,good]- rho[j,k]*oz[j,k,good]*v[j,k,good])/float((size(good))[1])
    qwu[k,j]=total(rwu[j,k,good]- rho[j,k]*oz[j,k,good]*u[j,k,good])/float((size(good))[1])
    qvu[k,j]=total(rvu[j,k,good]- rho[j,k]*oy[j,k,good]*u[j,k,good])/float((size(good))[1])
  endfor
endfor
nlev=64
qwvn=qwv/rurms2
qwun=qwu/rurms2
qvun=qvu/rurms2

save,filename='stresses.sav',qwvn,qwun,qvun
;
omega_0=1.e9/(28.*24.*3600.)                 ;28 days
omega[*,1:m-2]=1e9*(uu[*,1:m-2]/(2.*!pi*rsinth[*,1:m-2]))
omega=omega+omega_0

fo='(F9.6,2x,F9.6,2x,F10.5,2x,8(E11.4,2X))'
;openw,1,'T28.dat'
;printf,1,'     r    ', '     th   ', '     omega    ', '    u     ',$
;         '       v    ', '    quu  ', '   qvv  ', '   qww   ', '  qwv   ',$ 
;         '    qwu   ', '    qvu  '
;for k=0, l-1 do begin
;  for j=0,m-1 do begin
;    printf,1,r[k],theta[j],omega[k,j],uu[k,j],vv[k,j], $
;             quu[k,j],qvv[k,j],qww[k,j],qwv[k,j],qwu[k,j],qvu[k,j], fo=fo
;  endfor
;endfor
;close,1


thg=where(th GE -89 AND th LT 89)

xlength=0.26
xinterval=0.06
xpos1=xinterval
xpos2=xinterval+xlength

lev = 2.*max(abs(qvun))*indgen(nlev)/float(nlev-1)-max(abs(qvun))
print,'levels R_\phi\theta = ',minmax(lev)
cpos=[xpos1,0.15,xpos2,0.96]
bpos=[xpos1,0.055,xpos2,0.07]
contour,qvun[*,thg],x[*,thg],y[*,thg],/fi,nl=64,/iso,levels=lev, $
	title='!8R!D!7uh!8!N!6',pos=cpos,charsize=2,ystyle=4,xstyle=4
plots,circle(0,0,0.98),thick=3
plots,circle(0,0,0.62),thick=3
plots,circle(0,0,0.71),thick=2,li=2
oplot,[0,0],[-0.98,-0.62],thick=3,li=0
oplot,[0,0],[0.62,0.98],thick=3,li=0
colorbar,range=[min(lev),max(lev)],pos=bpos,/horizontal,xtickformat='(F7.2)'$
        ,xticks=2,xtickv=[min(lev),0.,max(lev)],yaxis=0,char=1.5
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength

lev = 2.*max(abs(qwun))*indgen(nlev)/float(nlev-1)-max(abs(qwun))
print,'levels R_\phi\r flux = ',minmax(lev)
cpos=[xpos1,0.15,xpos2,0.96]
bpos=[xpos1,0.055,xpos2,0.07]
contour,qwun[*,thg],x[*,thg],y[*,thg],/fi,nl=64,/iso,levels=lev,$
	title='!8R!D!7u!8r!N!6',pos=cpos,charsize=2,ystyle=4,xstyle=4
plots,circle(0,0,0.98),thick=3
plots,circle(0,0,0.62),thick=3
plots,circle(0,0,0.71),thick=2,li=2
oplot,[0,0],[-0.98,-0.62],thick=2,li=0
oplot,[0,0],[0.62,0.98],thick=2,li=0
colorbar,range=[min(lev),max(lev)],/horizontal,pos=bpos, $
         xtickformat='(F7.2)',xticks=2,xtickv=[min(lev),0.,max(lev)],xaxis=0,char=1.5
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
lev = 2.*max(abs(qwvn))*indgen(nlev)/float(nlev-1)-max(abs(qwvn))
print,'Total R_\theta\r = ',minmax(lev)
cpos=[xpos1,0.15,xpos2,0.96]
bpos=[xpos1,0.055,xpos2,0.07]
contour,qwvn[*,thg],x[*,thg],y[*,thg],/fi,nl=64,/iso,levels=lev,$
	title='!8R!Dr!7h!8!N!6',pos=cpos,charsize=2,ystyle=4,xstyle=4
plots,circle(0,0,0.98),thick=3
plots,circle(0,0,0.62),thick=3
plots,circle(0,0,0.71),thick=2,li=2
oplot,[0,0],[-0.98,-0.62],thick=2,li=0
oplot,[0,0],[0.62,0.98],thick=2,li=0
colorbar,range=[min(lev),max(lev)],pos=bpos,/horizontal,$
         xtickformat='(F7.2)',xticks=2,xtickv=[min(lev),0.,max(lev)],xaxis=0,char=1.5
;
;partvelvec,fx[*,thg],fy[*,thg],x[*,thg],y[*,thg],/over,fraction=0.1,length=0.08

!p.multi=0
end
